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Abstract 

The effect of interstitial hydrogen on the elastic properties of a-iron is investigated using ah 



Oh initio density functional theory. We find that while the overall strength properties are reduced by 

^ H, the effects are mainly due to the resulting local volume expansion. We use these concentration- 

^ dependent elastic moduli to model the effects of homogeneous strain fields on H solubility and 

make extrapolations to H-solubility in the spatially-varying strain fields of realistic dislocations. 

O 

^ We find that H is strongly trapped by dislocations of both the edge and screw variety, leading to a 

^ remarkable increase in the local H concentration near the dislocation cores. The strain dependence 

of the solution energy leading to accumulation of H near dislocation cores, as well as the reduction 
g of elastic moduli which our calculations predict, is consistent with numerous experimental studies 

which indicate that the trapping of H near dislocations plays a significant role in the embrittlement 
O of iron and steels. 

> 

00 
!>• 

cn 

o\ 
o 
o 



X 
J3 



1 



I. INTRODUCTION 



Hydrogen degrades the performance of many metals by lowering the failure stress, lead- 
ing to fracture at unpredictable loading condition^^. Several mechanisms have been pro- 
posed to explain the H-embrittlement of iron: the hydrogen-enhanced decohesion (HEDE) 
mechanisnpEI, H-vacancy effectd^, hydrogen-enhanced localised plasticity (HELP)'^'^^, or a 
combination thereof. However, none of these explanations alone seems to provide a complete 
explanation of the experimentally-observed features of H-embrittlement. 

In the HEDE mechanism, H weakens the cohesive bonds between the metal atoms, leading 
to failure at interfaces, as H exhibits increased solubility in the tensile strain field of a crack 
opening. Vacancies containing H can order themselves along critical slip directions, leading 
to fractur^^. Within the HELP mechanism, the onset of plasticity with loading occurs at 
a lower stress as a result of the H shielding of repulsive interactions between dislocation^. 
The increased H concentration at dislocation^ effectively reduces the dislocation-dislocation 
distances. Together with the experimentally-observed increase of dislocation mobilitj^^, this 
may explain the phenomenon of dislocation coalescence, and ultimately, crack advancement 
at reduced loads. The importance of increased H-concentration near dislocations and other 
low-energy trap sites was also indicated by in a recent experimental study of intergranular 
failure in steeP^. In that study, the fracture mode changed from ductile to brittle, as the 
amount of H located at these low-energy trap sites increased, indicating that a study of H 
at dislocations is important. 

The goal of this study is to achieve a microscopic understanding of the changes in the 
mechanical properties due to H, and their connection with the proposed theories of H- 
embrittlement. Focussing on a-Fe, we study the influence of interstitial H on the elastic 
properties, and the behaviour of H within strain fields representative of dislocations. To 
achieve this, we use density-functional theory (DFT) calculations to calculate the changes 
in the elastic properties of a- (bcc) iron, as a function of H-concentration and employ these 
results in a continuum elasticity theory description of the H-solubility in strain fields. In Sec. 
HI] we describe the details of our DFT calculations while Sec. IIIII describes the calculations 



of the modification of the elastic parameters as a function of H-concentration. Sec. IV 



discusses the solubility of H in strain fields, and we conclude with a discussion in Sec. |V| 
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II. DESCRIPTION OF FIRST-PRINCIPLES CALCULATIONS 



Our first-principles density functional theory calculations were performed using the 
Y^gpHSHIS] qq({q ^J^Q used the projector augmented-wave methocP^^^ and employed pseu- 
dopotentials considering the 3p electrons of Fe as valence electrons. The generalized-gradient 
approximation (GGA) in the PW91 parametrisatiorP^ was used for the exchange-correlation 
functional. All of our calculations were spin-polarised. Some of our calculations were re- 
peated using the Vosko-Wilk-Nusair^^ (VWN) spin interpoilation in PW91, or the Perdew- 
Burke-ErnzerhoP^ (PBE) exchange-correlation functional and are discussed in Sec. HID 



A plane-wave basis with a cutoff of 500 eV and supercell geometry were used. A F-centred 
k-point grid equivalent to 18x18x18 for the two-atom basis bcc unit cell was used for the 
Brillouin-zone sampling except for the 128-atom Fe cell where it was 20x20x20. Our re- 
sults for the lattice parameter (2.832A), bulk modulus (194 GPa), and magnetic moment 
(2.17 fis) are in good agreement with other first-principles calculations e.g. Refs. |2T|22] 
and with experiment I^HP^ The plane-wave basis and k-point sampling density used are 
more than adequate for converging the equilibrium properties, but necessary for converging 
elastic moduli to within less than one percent. The ions were relaxed such that the maxi- 
mum component of the force on each ion was less than 0.01 eV/Awhile the total energies 
were converged to within 0.01 meV. The lattice parameter for systems constrained to be 
cubic (see following section for a discussion of this approximation) was determined from 
fitting the total energies to the Murnaghan equation of statd^. At times, the correct lattice 
parameters for a tetragonal unit cell were desired; these were obtained from a quadratic fit 
over a two-dimensional grid of total energies. 

In this study, we focused our attention on tetrahedral interstitial H, for two reasons: one 
is that at zero stress we found the tetrahedral interstitial position to be 0.13 eV more stable 
than the octahedral position, in agreement with previous DFT studied^^*^, and secondly, as 
there are twice as many tetrahedral as octahedral sites per Fe atom (six versus three), we 
expect the former to dominate the mechanical properties at ambient temperatures. 
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III. ELASTIC PROPERTIES OF BCC-IRON 



A. Calculation of elastic parameters 



The total energy of a solid at zero stress and equilibrium volume Vq can be expanded 
about small strains using the Voigt notation (m running from 1 to 6): 

E{e) = E{0) + C^i^^^i + E C^Jkeie,ek + . . . . (1) 

ij ' ijk 

In this section we restrict ourselves to linear elastic behaviour, i. e. stress linear in the strain 
the first two terms of Eq. [TJ For the calculation of the elastic parameters in this work, cubic 
symmetry was enforced. This resulted in an error of less than 2%, as will be described at 
the end of this section. The number of independent elastic parameters Cij depends on the 
crystallographic symmetry. For cubic systems, there are three unique elastic parameters, 
Cii, Ci2, and C44, that describe the change in energy due to applied strain. The parameters 
Cij can be determined by applying suitable strain tensors and taking the second derivative 
of Eq. [1] with respect to the amount of strain. The combination -B = | (C^ii + 2C12) is the 
bulk modulus, found by applying hydrostatic strains and then doing a fit to the Murnaghan 
equation of state, whereupon the pressure derivative B'{P) is also found. The combination 
C = Cn — C12, and the coefficient C44 are found by applying volume-conserving strains. 
In the case of C, an orthorhombic strain 



SO 
0-6 
\0 6^(1 -6') J 

is applied, giving a total energy expression 
while, for C44, the strain is monoclinic: 



(2) 



(3) 



with a total energy 



Er, 



6/2 
5/2 
57(4-^2) J 

1 



(4) 



(S) = E{0) + ^C4^VS^ + O[S'] 



(5) 



We fitted our total energies for systems subject to nine values of strains S between ±3% 
to Eqs. |3] and [S] in order to obtain C and C44. To estimate the accuracy of the elastic- 
parameter calculations, extensive convergence tests were done. The quantities varied were 
the k-point density, cutoff energy, and percent strain (value of 6). The lattice parameter 
was fixed at 2.832 Afor all tests. An example of these convergence tests for C44 is shown 
m Figjlji, for 6=3% (the convergence tests for C are similar). The values of cutoff energy 
and k-point density of the circled point (500 eV, 18x18x18) were used in (b) and (c) to 
check that the applied strains are within the linear regime. As a result of our extensive 
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FIG. 1: (a) Convergence of C44 as a function of k-point density and energy cutoff for the two-atom 
bcc-Fe cell for a range of \5\ < 3%. The error bars are deduced from the least-squares fit to Eq. [5| 
The circled point corresponds to the calculational parameters used in this study, (b) C44 and (c) 
C fitted between different ranges of ±6. 

convergence tests, we observed less than a 1% variation in the elastic parameters upon going 
from 0.5-3% change in S. For the calculations presented in the remainder of this study, we 
chose a value of 3% for 6 in order to avoid numerical instabilities and/or the need for very 
high presicion total energy evaluations. 



B. Elastic parameters of FeH 



The volume expansion of the crystal lattice by H in the tetrahedral interstitial site in- 
troduces small tetragonal distortions, breaking the cubic symmetry of the Fe host lattice. 
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An H-orientation such as the one shown in Fig. |2] causes an expansion of the Fe nearest- 
neighbours radially-outward from H (along the red lines) This expansion, projected onto the 
cube axes, is of an identical amount in the a and c direction, and greater than the expansion 
along b. The tetragonal distortion, with lattice parameters a = c ^ b, increases the number 




a 



FIG. 2: Location of a H-interstitial on the face between two bcc unit cells. Lines towards nearest-Fe 
neighbours (blue circles) are shown. 

of unique elastic parameters from three to six: C12 = C23 7^ C13, Cu = C33 7^ C22, and 
C44 = 6*66 7^ C55. The variation among the no- longer equivalent elastic parameters depends 
on the degree of tetragonal distortion, generally becoming more pronounced with increasing 
H-concentration, unless there are two or more H in the supercell, in which case the amount 
of distortion depends also on their relative positions. Generally, H would be expected to be 
randomly distributed, resulting an isotropic lattice distortion. As well, considering a cubic 
unit cell facilitates the scale-bridging with mesoscopic approaches such as e.g. finite-element 
schemes. 

The error due to imposing cubic symmetry can be estimated by explicitly calculating 
the (formally-degenerate in cubic symmetry) elastic parameters, for tetragonal symmetry 
for a configuration with lattice parameters a = c 7^ 6 by straining the system in all three 
cartesian orientations, all the while maintaining a cubic unit cell. The reason that the results 
are not equivalent in all directions even under the assumption of cubic symmetry is because 
of the local internal tetragonal symmetry around the H atom. Nevertheless, the variation 
amongst the three orientations was found to be small: in two of the orientations the elastic 
parameters were identical, as expected, while those from the third differed by about 4%, 
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making the standard error less than 2% for the highest concentration (11%) studied, where 
the tetragonal distortion 1 — h/a was -0.004, where the two H were oriented such that the 
distortion is double that caused by a single H. For all of our data above and including a 
H-concentration of 4%, the formally-degenerate elastic parameters were explicitly calculated 
by straining in the three cartesian directions and the results presented are the average values 
along with the associated standard errors. For lower concentrations, the spread amongst the 
three directions was < 1% and corresponding errors presented are from the least-squares 
fitting of the total energy as a function of strain. 

In addition, for the system with 1 — h/a = —0.004 (at 11 at.% H), the elastic param- 
eters for an explicitly tetragonal lattice were calculated using the strain and total energy 
expressions given in Ref. EHl The differences from the simplified case of the cubic unit cell 
with distortions applied in all three orientations were small, less than 2% for the largest 
H-concentration studied. Table [T] summarises the average values of the cubic and tetragonal 
elastic parameters along with the associated standard errors computed from the spread in 
the values for the different orientations. Only C12 gave a large spread for the tetragonal 
cell (separately, the values were C12 = 6*23=1.38, Ci3=1.64). That the agreement is good 
is not surprising, since the volume of the cubic cell was made to be the same as that of 
the tetragonal cell and the variation in the elastic parameters arising from different lattice 



parameters (Sec. IIIC) in the tetragonal cell mostly averages out in the cubic cell. 





cubic 


tetragonal 


Cn = g {Cu + C22 + C33) 

C12 = g (C12 + Ci3 + C23) 

C44 = 1 (C44 + C55 + Cee) 


245±2 
142±2 
92±2 


239±6 
147±15 
93±2 



TABLE I: Elastic parameters (GPa) for the highest-considered H-concentration (11 at.%) for a 
cubic or tetragonal (l-b/a=-0.4%) unit cell, both having the same volume. For the cubic cell, 
the strains were still performed in the three, formally-degenerate directions, and averaged, as the 
results were different due to the lower symmetry locally around the defect. 
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C. Non-linear elastic behaviour of pure Fe 



In order to achieve a comprehensive description of the influence of volume-expansion due 
to interstitial H, we also investigated the nonlinear elastic behaviour of the elastic parameters 
as a function of hydrostatic strain. Equivalently, this is the effect of a volume expansion on 
the linear elastic parameters Cij. The modified second-order elastic parameters are given by 
Birch^^ in terms of the second and third-order elastic parameters Cij and Cijk to linear order 
in the applied hydrostatic strain. Instead of calculating the third-order elastic parameters, 
we directly determined the variation in the second-order elastic parameters as a function of 
small applied hydrostatic strain such that the assumption of linear elasticity remains valid. 

Following Wallace [30], the total energy at a volume V produced by strains applied at 
a reference volume V^'^'^ away from equilibrium, is modified from Eq. [l] to include a term 
first-order in strain, corresponding to hydrostatic stress aij: 

E{V, e) = E{V'f, e = 0) + V''^ ^ a,,e,, + C^jki{V'-'^)eijeM. (6) 

ij ijkl 

In Eq. |6| (as in Eq. 2.37 in ^30j) the strains are infinitesimal strains evaluated with respect 
to V^"^-^ . Because the stress is hydrostatic, it can be written as aij = —P6ij. Expanding Eq. 
[6]and allowing the strains represented by Eq. [2] and |4] to be applied with respect to V^'^-^ , 
we obtain the expressions, modified from those of Eq. |3] and [5} 

Eortho.iP, S) = E{P, 0) + (C - P) ^ Q[^4]^ (7) 

EmonoXP, 5) = E{P, 0) + ^ (Cu - V'^S^ + 0[S% (8) 
At each value of 77 ranging between ±0.05, the value of 6 was varied between ±0.03 and a fit 



was made to second-order in 6, as in Sec. Ill A The elastic parameters C and C44 are again 



extracted from the coefficients of 6"^. However, in order to separate Cu and C12 from C, the 
bulk modulus needs to be adjusted to take into account the nonzero pressure, by an amount 
equal to its derivative found from the fit to the Murnaghari^ equation of state (we found 
B'{P) = 5.42), multiplied by the pressure at the corresponding value of hydrostatic strain 
rj. The variation of the elastic parameters with hydrostatic strain follows a linear trend, as 
can be seen from Fig|3] A linear least-squares fit to the small-strain region (|?7| < 0.02) of 
these curves resulted in the following dependence of the elastic parameters on hydrostatic 
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strain: 



B = 


(194- 


lOeir]) GPa 


Cu = 


(284- 


1640r/) GPa 


Cl2 = 


(149- 


798ri) GPa 


Cm = 


(105- 


771ri) GPa 



Previously-reported first-principles calculations of the elastic parameters of bcc-F^^^^ 
show a spread of approximately 10%. Some of these results, including dependences on hy- 
drostatic strain, if reported, along with values from an experiment performed at 4 Kelvin, 
are shown in Table |TT] along with our results. We did a fit of the others' data to obtain the 
reference values listed in the table. The experimental data for the hydrostatic-strain depen- 
dences of the elastic parameters is sparse. The few experimental data points existing in the 
literature are superimposed onto the calculated results by Sha and CoherP^ in their paper, 
and these agree well with their calculations. Our derivatives do not deviate appreciably 
compared to the others' results. 



method 


Cii(GPa) 


Ci2(GPa) 


C44(GPa) 


PAW-GGA (present) 


284(-1640) 


149(-798) 


105(-771) 


LMTO-GGAED 


303(-1282) 


150(-813) 


126(-604) 


PAW-GGAEH 


271(-1228) 


145(-535) 


101(-454) 


PP-GGAISI 


289 


118 


115 


FP-LAPW^ 


285 


139 


100 


expt.ISS 


245 


139 


122 


expt.l^ 


240 


136 


121 



TABLE II: Comparison of elastic parameters (and their hydrostatic strain rj dependence in paren- 
theses, if available) obtained in the present study with other calculations and experimental results. 

D. Dependence of elastic parameters on H-concentration 

The elastic parameters of pure Fe serve as a starting point for determining the influence of 
interstitial H atoms on the elastic properties. We implement the variation of H-concentration 
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-0.05 -0.025 0.025 

hydrostatic strain rj 



0.05 



FIG. 3: Calculated elastic parameters for pure Fe as a function of applied hydrostatic strain 
rj = AV/Vq where Vq is the equilibrium volume of bcc-Fe. 

in the supercell approach of our calculations by (i) increasing the size of a supercell containing 
one H atom, or by (ii) adding a second H atom to the same supercell. In addition, we changed 
the symmetry of the supercell or the position of the second H atom relative to the first in 
order to alter the ordering of the H atoms within the Fe host lattice. The sizes of supercells 
and number of H atoms used to achieve various H concentrations are listed in Table IIIII 
Figure |4] shows the variation with concentration of the elastic parameters. Additional data 



points at the same concentration correspond to differently-ordered structures (see Table III). 
There is a clear trend of decreasing elastic parameters with concentration. We fit a line to 



each set of data, which we will make use of in Sec. IV This line was made to include the 



zero-concentration value and excluded the outliers at 0.8% and 10%. We ascribe the outlying 
values of Cu and for the 10% concentration to the fact that this supercell had the largest 
stress associated with it. This was the case because of the relative orientation with which the 
two H atoms were set up, and also to the proximity to the periodic H images in the smallest 
dimension. The C44 at 0.8% (and to a lesser degree at 2.7%) lies outside the line of best fit 
and warrants further investigation as to whether it is a function of the precision used in the 
calculation or a genuine effect. We also repeated some of these calculations using the VWN 
corrections, and also with the PBE functional. They resulted in slightly different equilibrium 
lattice parameters (order of 0.3%), and correspondingly, different absolute moduli, but the 
slope of the Cij versus concentration remained the same. 

Calculations of the modification of elastic parameters of Fe by H has been reported 
recently in Ref. [36l Their underlying lattice was not perfect Fe but a crack with tens of 
thousands of atoms and an atomistic model was used for the calculations. The data covered 
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at. % H 


n(H) 


kxlxm 








1x1x1 


0.8 


1 


4x4x4 


1.8 


1 


3x3x3 


2.7 


1 


3x3x2 


3.6 


2 


3x3x3 


4.0 


1 


2x2x3 


5.3 


1 


3x3x1 


5.9 


1 


2x2x2 


5.9 


1 


2x4x1 


7.7 


1 


1x2x3 


10.0 


2 


3x3x1 


11.1 


2 


2x2x2 



TABLE III: H-concentrations and how they were achieved based upon the number n(H) of H and 
supercell dimensions in multiples k, I, and m of the lattice parameters of the two-atom bcc Fe 
cell. F the highest concentration (11.1%), three configurations, with H-H spacings within the same 
supercell of 4.98, 3.66, and 2.31 Awere calculated. 

up to 6 at. % H and while the elastic parameters decreased with H initially, the effect levelled 
off with increasing H and generally was far weaker than here. 

Figure [5] shows that the volume of the supercell changed linearly with H concentration, 
and we fit a line through this data. These volume changes may be related to the strain 
parameter t] of FigjSj and Eqs. |9] to assign the changes in the elastic parameters as arising 
from volumetric effects or other effects. 

The volume-induced changes in elastic properties are shown, together with the total 
effects orignally displayed in Fig. |4], in Figure |6j 

Single-crystal Fe samples with H are difficult to prepare while the measurements of the 
Cij do not directly relate to the strength properties of the material. Most often, the samples 
are polycrystalline, and the usual measurements performed are those directly related to 
stiffness, tensile strength, and hardness, in other words the bulk. Young's and shear moduli. 
In addition, for higher-length-scale simulations such as finite-element calculations of realistic 
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FIG. 4: Calculated elastic parameters for bcc Fe as a function of H concentration. 
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FIG. 5: Relative volume expansion as a function of H concentration where Ay is the volume 
difference of the system with H compared with the pure Fe volume V'^ . 



materials, the polycrystalline quantities are more practical than the single- crystal ones. 
Motivated by this, we used combinations of the elastic parameters Cij, to derive various 
elastic moduli, describing different types of stress-strain responses. Shown in Fig. [7] are the 
bulk modulus, and the polycrystalline averages of the Young's and shear moduli, where the 
average was performed according to the expression by Hill, which is an average of the Voigt 
and Reuss bound^. As in Fig. [6| the volumetric-strain-induced changes in the moduli for 
pure Fe are also shown. 

The literature is extremely scarce on relevant measurements of elastic moduli of Fe with 
H. A rare exception is Ref. [3B1 where the effect of H on the shear modulus of polycrystalline 
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FIG. 6: Variation of elastic parameters from dft and from volume dependence of pure iron's elastic 
parameters 

Fe was studied. They extrapolate their results to 1 at.% H and predict a decrease of 8% in 
the shear modulus at a temperature of 100 K, in agreement with our calculations. 

By taking the difference of the two curves in Fig. |6] or Fig. [7], we remove the contribu- 
tion of strain of pure Fe from the elastic parameters calculated at the equilibrium volume 
corresponding to each concentration of H. The resulting graphs, in Fig. [8} show only the 
electronic effect of H at each concentration (and implicit corresponding volume). These 
graphs show a clear trend in the elastic regime towards a stiffening, tensile strengthening, 
and hardening of the system with increasing concentration of H. Despite the overall soft- 
ening of the material with increasing concentration, it is clear, that the effect of H leads 
to less softening than would be expected from the volume expansion alone. The separation 
of solute effects on the matrix elastic moduli has been also considered in Ref. [39] for H 
in Nb. Instead of parametrising the volume effect by the second term of their Eq. 4, we 
explicitly calculated it. Otherwise it's a straight line passing through our volume-derived 
elastic moduli curves. Similar analysis was done in Ref. SOl 
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FIG. 7: Variation of polycrystalline moduli from DFT data (Fig. |4]), and from the volume depen- 
dence of pure iron's elastic parameters 




4 6 
) atomic H concentration 



FIG. 8: Removal of volume-dependence of elastic moduli (subtracting DFT data from V- 
dependence in Figs. [oItI). 



IV. STRAIN-DEPENDENCE OF H SOLUBILITY 



In this section, we describe the effect of different types of strain on H solubility. Our goal 
is to examine the effect of H on typical dislocation strain fields and conversely, the latter's 



14 



impact on H solubility. In order to avoid the computationally very demanding explicit 
treatment of dislocations in correspondingly large supercells, we have instead examined the 
solubility of H in different types of homogeneous strain fields which we later used to estimate 
the H-solubility in realistic dislocation elastic strain fields. 

Specifically, we have calculated the heats of solution for tetrahedral H in supercells subject 
to hydrostatic, or uniaxial or shear strain. These applied strain tensors mimic the strain 
fields of different dislocations as discussed in detail below. 

To arrive at the heat of solution, W, at each value of applied strain, we subtracted from 
the total energy of the cell with N Fe atoms and n H atoms, EpgH„, the total energy of 
a pure iron cell with N Fe atoms, and half the energy of a free H2 molecule multiplied by n: 

Win) = EpeH^ - Epe - u^Eh, . (10) 
In all calculations, the cells with H were constrained to be cubic in their unstrained state, 



which is a good approximation, as described earlier, in Sec. Ill B 



Equation [10} however, can be interpreted in different ways depending on the strain or 
stress state of the pure Fe reference system with respect to the system with H. Therefore, 
as the definition of solution energy is ambiguous, we have defined four regimes on which 
to concentrate as the strain of the FeH system is varied: (i) one where the pure Fe volume 
matches that of the FeH system at each strain, (ii) where the two systems have the same 
stress tensor at each strain, and (iii) where they have the same strain tensors. 

A. Parametrisation of the H-solution energy 



Using our DFT-obtained material parameters (Table IV), we have parametrised the so- 



lution energies (Eq. 10) as functions of strain and H-concentration and we find them to 
accurately describe our DFT results, to a large extent quantitatively. We label all quantities 
with an a; or superscript for the FeH or pure Fe system respectively to denote which system 
they correspond to, for example, for the equilibrium volume of the FeH system and 
for that of pure Fe. For a system Y ((=0 for Fe, x for FeH), the total energy at an arbitrary 
reference volume V^'^'^ and infinitesimal strain tensor e with respect to V^'^'^ is (c./. Eq. [6| 

{V^-f, e) = E'^iV^-f, 0) + + ^ E CrMV''^)^^je,,, (11) 

ij ijkl 
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FIG. 9: DFT-calculated solution energies along with parametrised form for (a) hydrostatic strain, 
(b) uniaxial strain, (c) shear strain, for different H-concentrations. 

where the first term is the total energy at the reference volume V^'^-^ , the second term is the 
stress part which is non-zero if V^'^-^ is not the equilibrium volume V"^ of Y, and the last 
term is the elastic deformation energy about V^''^-^ . Explicitly, the first term is equal to 



1 ...... ..Y^2^^E^ 



E^iV'^f^ 0) = 0) + - - r^) 



V=VY 



E^{V^,0) + ^{V-f-V^)\ (12) 
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where -B^ is the bulk modulus of Y at . The second term, assuming a hydrostatic stress 
corresponding to a pressure to achieve the reference volume is 



Y 



QYyref 



(13) 



where the relation = —B^ (V^^^ — V^) /V^ has been used. This relation for assumes 



linear elasticity holds (5^ doesn't vary) but the contribution of Eq. 13 is of 0{x^) and as 



such this assumption is reasonable for low concentrations. In the third term of Eq. 11, the 
elastic parameters are evaluated at the volume V^'^^ and as such, they are equal to 



dV 



1 F)r'^ 



(14) 



dr] 

where the volume derivative has been rewritten in terms of the hydrostatic strain rj = 
^yref _ yY^ /V^ . Despite their generality, Eqs. 



12 



14 



are difficult to work with for the FeH 



system, as we do not have explicit data on the strain-dependence of the C^ji^i- Therefore, 
when dealing with the FeH system, we take the reference volume as V^'^-^ = V^. The volume 
may be expanded about V^, for small x as = V'^ + xdV/dx. 

We now list the expressions pertaining to FeH. These will remain the same throughout 



for all applications considered. The first term of Eq. 11 [i.e. Eq. 12) is equal to 



E^(l/^,0) = E^{V^,^)+x 



1 



0\2 



+ 2 (^'' - 
dE 



dE 
dx 
d^E 



+ x{V'' - V^) 



d^E 



dxdV 



dV 



V°,x=0 

dx dx dxdV 



1 2 d^E 

^2^ 



+ 



B 







2 



dV\ 



X- 



1 2^ 

+ 2^ dx^-' 



(15) 



2^0 \ dx J 

where here and in what follows, the limits of evaluation of the derivatives are omitted for 



clarity. In Eq. [15} the first derivative with respect to volume does not appear because it 
is equal to zero as the evaluation limits pertain to pure Fe. The second term of Eq. 11 



i.e. Eq. 13) is zero because is the system's equilibrium volume. For the third term, the 
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elastic parameters (Eq. [14]) are 



^ijkl = ^ijkli^^) = ^ijkl "I ^^^^ 



where the derivative represents the variation not only with respect to concentration about 
but with strain implicitly included as well. These derivatives are the slopes of the lines 
of best fit from Fig. |4} 

After defining the relative volumes, stresses, and/or strains, of the FeH and Fe systems 



corresponding to each of the cases outlined above, we use Eqs. 11 and 10 to derive strain- 
dependent solution energies for different relative strain and/or volume conditions for the Fe 
and FeH systems. 



1. same volumes 



The case of equal volumes would refiect the energy to be gained or lost by inserting \H2 
into a strained Fe lattice without further changing the average volume per Fe atom. Since 
the reference volume of the FeH system is that is also the reference volume of the pure Fe 



system. This results in a contribution from Eq. 13 We choose the strains to be e = e^, the 
FeH system strains. The pure Fe system, at a reference volume V^, applying the equations 



of the previous section (Eq. 11 and following), is described by a total energy 

2 



E 







2^0 



^ ' i ijkl 



(17) 



where the elastic parameters are, applying Eq. [T^ 



1 



1 9^0,, dV 
-X- 



drj dx 



(18) 



The solution energy, obtained by collecting all the contributions to Eq. 11 and then using 
Eq. 10, neglecting all terms of order except for those appearing in Eq. 15, because we 



want to explain only the discrepancies at zero strain between the various concentrations, is 
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Quantity 


value 


A E 


0.22 eV 




194.2 GPa 


A V 


4.5 A3 




11.352 A3 


clB^/drj 


-1061 GPa 


dC?^ /dT] 


-1508 GPa 


dC%/drj 

ILI 1 


-885 GPa 


dClJdri 


-677 GPa 


dB^^^dx 


-153 GPa 


dC{f/dx 


-316 GPa 


dC\^/dx 


-72 GPa 


dClf/dx 


-141 GPa 



TABLE IV: Parameters used in the expressions 19 - 23 Also used in subsequent sections. All 
derivatives are evaluated at and x = 0. 



found for a cubic system to be equal to 

W/n 



+ 



2N dx^ 



+ 2 
+ 4 





dx 


drj 




dCif 


dC'u 


N 


dx 


drj 






dC^^ 



dxdV 

I \ XX 



+ 2 + 2) 



'ZZ 



yy ' ^xx^zz 



' ^yy^zz; 



(19) 



dx dr] 

In the above, we have used the definition of concentration as the ratio of the number of H 
atoms, n, to the number of Fe atoms A^: x = n/N, to rewrite the concentration derivative of 
V in terms of the volume expansion AV per H and A^: ^ = 4^4^ = AVN and similarly for 

dx an ax •' 

the concentration derivative of E, ^ = AEN. The ratio V'^/N is the volume per Fe atom at 
equilibrium, a constant. Numerical values of all quantities appearing in the parametrisation 
are listed in Table 



IV 



A£'=0.22 eV was obtained by calculating the solution energy 
in the largest unit cell (128 atoms) studied. The parametrised curves in all of the figures 
were generated by assuming a vanishing concentration, such that the contribution from the 
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x{%) 


A^+-^0(eV) 

Ziv ax^ ^ ' 




sum 


DFT T^(x)L x_n(eV) 


0.8 


0.246 


-0.023 


0.223 


0.229 


1.8 


0.238 


-0.048 


0.190 


0.192 


2.7 


0.249 


-0.066 


0.183 


0.190 


3.6 


0.251 


-0.089 


0.162 


0.175 


4.0 


0.250 


-0.093 


0.157 


0.173 


5.3 


0.278 


-0.123 


0.155 


0.177 


5.9 


0.255 


-0.149 


0.106 


0.136 



TABLE V: Energy terms from Eq. 12 contributing to the intercept at zero strain for the case of 



equal volumes (Fig. [9j). The first two energy columns are from expression 15 and were calculated 
using DFT as described in the text. Their sum is under the heading sum. The last column is the 
DFT data. 



x-dependent parts of Eqs. 091- 23 



IS zero. 

For the case of hydrostatic strain e = e^x 



AE H — — + xAV- 



" ^yy 

3B^AVe + 



zz, Eq. 19 becomes 
9 fV^dB*"* 



2N dx^ ' ' dxdV 
For uniaxial strain e = e^x and shear strain e = e 



2 \ N dx 



dB^ 
di] 



AV ' "2 



(20) 



xy 



^yx 



the resulting expressions had 



the form 

y^uniax.l^ 



AE 



X d^E d^E .AT. 

+ xAV—^ - B^AVe 



2N (9x2 



dxdV 



1 ^V^dCif _ dCl^ 
2\Ndx drj 



AV 1 



mishear I ^ = + + xAV^, 



(21) 
(22) 



2N dx^ dxdV ' \ N dx drj 

Having a combination of hydrostatic/uniaxial/shear strain would involve combining expres- 
sions 



20 22 as can be seen from Eq. 19 in any case. 
We performed DFT simulations with the Fe and FeH systems at the same volume and 
varied the amount and type of strain. We subtracted the two sets of energies and the energy 
of ^H2, as in Eq. 



10 



20 



22), using the DFT-calculated 



The parametrised forms (Eqs. 
parameters listed in Table |IV[ accurately describe the DFT-obtained solution energies, as 
can be seen from Figure |9} Minor variations are seen for the calculated solution energies 
at zero strain among the various concentrations, and this can be accounted for entirely by 
evaluating the terms x^ + x^^-^-Mj + \x'^'^ (c-/ Eq. 15 and the normalised form in 
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Eq. 19) explicitly at each concentration studied. The cross-term ^^^y was calculated at the 



point x=1.8% and this value (-0.557 eV/A) was used throughout. Using a graph similar to 

Fig. |9^, but with the horizontal axis in terms of the strain e° of pure Fe, + 

was predicted using the zero-strain data (equal to a volume of V^) from there. This can be 



V 



The sum of these 



done since W{x)\^o^q = AE + 2^^. These are given in the Table 
two terms, stated in the third column, agrees very well with the DFT data given in the last 
column. The cross-term was also calculated using data from the 5.9% concentration and it 
was a bit smaller (-0.426 eV/A) than the one used for Table [v[ leading to a contribution 
x2^^=-0.114 eV and thus a total E^(l^^', 0) =0.141 eV, at a; = 5.9%, thus accounting 
for the minor deviations at the higher concentrations between the two last columns of the 
table. The deviation of the intercept from the ideal value 0.22 eV for all the graphs shown 
in Fig. |9] can thus be explained. 

A derivation of the solution energy as a function of strain and concentration was pre- 
sented in Ref. HU In that work, the authors presented a systematic expansion of the total 
energy in terms of strain and concentration. They find, as we do, that the derivative of 
pressure with respect to concentration yields a contribution to the solution energy that is 
linear in hydrostatic strain. Using their expressions to express our concentration and strain- 
dependent elastic parameters requires considering up to the third derivatives of the total 
energy with respect to concentration and strain, which is the same as our concentration and 
strain-dependent elastic parameters. 

2. same stresses or same strains 



Another possibility of evaluating the solution energy Eq. [TO] is to compare the FeH and 
pure Fe reference systems with both having the same stress or strain tensors. We choose to 
therefore expand the total energy of each system about its equilibrium volume. 

For the two systems to have the same stresses (T„, we must have, in the Voigt notation, 

rr — \^ r'O ^0 — \^ r'^ 



k 



for m,n,k from 1 to 6. In the case of hydrostatic strain ef^ = Jjj-e, this leads to e° = e^o- In 



the case of the systems instead having equal strain tensors, we simply have e°- = 



^3^ 



As in the previous example, we apply Eq. [15] and we obtain, as a function of FeH 
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hydrostatic strain 
solution energy per H: 



Sije and shear strain ei^^^ 



5, the following expression for the 



+ 



9 



B''AV± 



N dx 



(23) 



The upper (lower) sign is for the case where the strains (stresses) are equal. Eq. 23 



IS 



sufficiently general such that expressions can be found for mixed cases, e.g. equal hydrostatic 



stresses but equal shear strains. The expression 23 is compared with DFT data in Figure 



10 The case 5 = and varying e for equal strains is shown in (a) while (b) has equal strains 
as well but here e = and 6 varies. In (c) we compare the case of equal stresses for e = 



and 5 varying. All the parameters were taken from Table IV, and converted to eV/A'^ prior 
to using. 

The higher curvature for the 2.7% concentration in (c) occurs because at this concentra- 
tion, the C44 data (Fig. |4]) lies below the line of best fit, thereby giving a larger contribution 
to the total derivative. 



As with the previous case of equal volumes, the intercept at zero strain in Fig. 10 



can 



be accounted for by explicitly calculating the zero-strain terms appearing in Eq. 23 The 
comparison, given in Table |VI| is very good for small concentrations and deviates for the 
larger one, by a similar amount as the case of constant volumes. As in that case, discussed at 



the end of Sec. IV A 1, the deviation is attributed here as well to variation of the cross-term 



with concentration. 



x(%) 


Table 


V 


sum 




sum 


DFT Ty(x)|,.^o(eV) 


2.7 


0.183 


0.030 


0.213 


0.220 


4.0 


0.157 


0.045 


0.202 


0.211 


5.9 


0.106 


0.068 


0.174 


0.202 



TABLE VI: Energy terms contributing to the intercept at zero strain (Fig. 10) for the case of equal 
strains or stresses. The first energy column is the same as that in Table [V] while there is now an 
additional term (second energy column in this table) appearing. The sum of the first two columns 
is under the third column, sum, which should be compared with the DFT data in the last column. 
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FIG. 10: DFT-calculated solution energies per H atom for various concentrations along with 
parametrised forms (Eq. [23]) for a reference pure Fe system with (a) same strains, hydrostatic 
only (b) same strains, shear only, and (c) same stresses, shear only. The strain on the horizontal 
axis is for the FeH system. 

B. H-solubilities in realistic dislocation strain fields 

Here we will look at the variation of H-solubility as a function of strain direction. As well, 
we aim to visualise the H-solubility in a typical dislocation strain field. We will focus on the 
particular case of a screw dislocation lying along (111), which is the predominant orientation 
in a bcc metaP^. The screw dislocation strain field was modelled using isotropic dislocation 
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theorj^. In cylindrical (r, 6, z) coordinates, its only nonzero component is = ^ze = 
where b = 2.45A is a Burger's vector equal to the length of ^(111). 

All of the figures and equations for shear strain (Figs. 10d-c, Eqs. 22 23) look the 
same whether the shear strain is e^y = e^^. = 6, e^^ = ef^ = 6, or e^^ = e^y = 6. For a screw 
dislocation line situated along (001), the nonzero components are e^^ = ef^, = e^^ = ef^ = 
Therefore, referring to Fig. [10} we see that H is trapped towards a screw dislocation core in 
the case where the pure Fe reference system has the same shear strain as the FeH system 



(Fig. 10 d). We aimed to see whether the repulsive behaviour of H by a screw disclocation 



at constant volumes (Fig. |9p) can be altered, or even reversed in orientations other than 
the coordinate system defined by the cube axes. 



The solubility equations [20p3] were examined to see whether they yield different values in 
rotated coordinates corresponding other orientations of the dislocation line. Depending on 
the choice of coordinate system, it is, in general, possible to get stronger or weaker variations 
for the solubility than those shown in Fig. |9} Already another combined atomistic/DFT 
study has shown that H-trapping in bcc Fe is very dependent to the direction of applied 
sheai'^. 

For the case of purely hydrostatic strain, we find that the solution energy was independent 



of the coordinate system chosen. This is because in our model (Eq. 19 ) both the term linear 
in strain as well as the quadratic term can be factored in terms of the trace e^x + ^yy + 
e^z- Since, as can be shown, the trace of the strain tensor is invariant under any unitary 
transformation, the solution energy for a hydrostatic distortion is thus rotationally-invariant. 
In the case of the screw dislocation, there is no such restriction; it is certainly possible to 
have different coefficients, some leading to H-trapping near a screw dislocation core. 



Figure 11 shows the variation of the solution energy for a fixed, 2%, shear strain where Eq. 



22 has been rotated onto other axes. By rotating the entire coordinate system, {x,y, z) — >■ 



{x' , y' , z'), other directions z' for the dislocation line (legend labels) are achieved. The label 
6 on the horizontal scale defines y' and hence the shear direction. As the directions (110), 
(112), and (111) are mutually perpendicular, the first two intersect the latter at ~ 35° 
and 6 = 7t/2 respectively. We expressed the strains ex'z' = ^z'x' in the primed system in 
terms of the unprimed system and inserted these values into the general expression for the 



solution energy Eq. 19 Note that only the e^,^, = e^,^, components are considered. For the 
full description of the screw dislocation, the e^,^, = e^, , components are also active. The 
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> 
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0.22 -■ 



0.215 - 



0.21 




FIG. 11: Variation of solution energy for 2% shear strain (Eq. 22) for a dislocation line/Burger's 
vector oriented along different directions (legend). The direction of the dislocation line/ Burger's 
vector is, black: [001], blue dashed: [111], red: [110], green dashed [112]. The horizontal axis is the 
variation of the solution energy as a function of the choice of x' and y' axes (and hence the shear) 
within the plane normal to the Burger's vector. 9 = has been arbitrarily chosen to correspond 
to the y' = {khO) direction. The dotted line is the unstrained value of the H solution energy. 

latter components result in the same curves as those depicted except with a 7r/2 phase shift 
such that the average value, corresponding to all components of a screw dislocation, is flat, 
independent of shear direction. In contrast with the study [S], we do not find asymmetry 
of our solution energy W about zero strain as a function of applied shear strain. To do 
otherwise would mean acquiring linear terms, which are not possible with our model. 

For each point on a grid in the plane perpendicular to the dislocation line, the solubility 



corresponding to the complete strain tensor was calculated from Eq. 19 in the case of 



equal volumes or analogous equations (not given due to their complexity) for the cases of 
equal stresses and equal strain. One caveat concerning the case of equal stresses is that we 
are not modelling the solubility of H in the stress field of a dislocation but its solubility 
at individual stresses at each point r, such that the material is not self-consistently in 
mechanical equilibrium. However, this is a first approximation in a self-consistent treatment 
of H-solubility SSI The results, shown in Fig. [T2| show that H is attracted to the core for 
equal shear strains or equal volumes and repelled for equal shear stresses. 



25 



> 



0.35 
0.3 

0.25 
0.2 

0.15 



I I 1 1 1 1 f 1 1 


' 1 ' 


- 


s. 


same stresses 




- \ 


same volumes 






■ — ■ ■ same strains 


- 














",1,1,1,1, 


1,1, 





3 



6 7 

r(A) 



10 



FIG. 12: Variation of H-solubility energy as a function of distance from the (111) screw dislocation 
core for the cases of equal stresses, equal volume, or equal strains for the FeH and pure Fe systems. 
The dotted line denotes the zero-strain solubility. 

C. Concentration of H around strain fields 



Motivated by the strong strain-solubility effects, we aim to calculate the local modification 
of the H concentration by strain fields. For the purposes of illustration, we consider only 
the case of equal volumes in this discussion. We consider the change in free energy F as a 
result of adding H to Fe under strain. The free energy is 



E 



FeH 



TS, 



(24) 



where T is the temperature, and 5* is the total entropy of the FeH system. In order to isolate 
the effect of adding H to the already-strained Fe lattice, we assume that the entropy of the 
H is independent of that of the Fe lattice: i.e. S = Sh + Spe- EpeH is the total energy of 
the (strained) FeH system. It can be rewritten in terms of the solution energy W {c.f. Eq. 
10 ) as EpeH = Epe + W — /in, where we have set the energy per H as — /i + |i?2, in contrast 
where it was equal to 



with Eq. 
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The free energy is minimized with respect to n: 





6E 
Sn 



dW _^dS _ 
dn dn ^ 



(25) 



We first focus on the case of a hydrostatic strain field. For simplicity we truncate Equation 
20 from the previous section to linear order in strain: 



W = W{e) = n{AE - 3B°AVe). 
26 



(26) 



The expression for the entropy can take on various forms, in the presence of various forms 
of H-H interactions, all equivalent at lower concentrations. In what follows, n is the number 
of H, Nq is the number of matrix (Fe) atoms, = QNq is the total number of tetrahedral 
interstitial lattice sites, and c is the occupation concentration defined as c = n/N. The sim- 
plest form for the entropy, 5*0, arises from an unrestricted concentration of H, corresponding 
not to a lattice but to an ideal gas. In this case, the entropy derivative is ^ = —k\n^. 
The concentration yielded by Eq. |25] for the case of hydrostatic strain up to linear order in 
the strain takes the Maxwell-Boltzmann form of 

c = coexp [3B^AVe/kT] 

where Cq = exp [{—AE + fi)/kT]. This reproduces the result from the Cottrell picture of 
enhanced concentration around a dislocatiorJ^l. A more realistic form of the entropy takes 
into account the presence of a lattice and enforces a maximum occupancy of one H per 
tetrahedral site. In this case the entropy is Si = —k [^In-^ + (A^ — n)ln (l — -^)] and the 
entropy derivative is ^ = — /cln-^y^. The concentration using this form follows Fermi-Dirac 
statistics: 

c = {exp [AE - 12- W^AVt/kT] + l}"^ 

However, the maximum occupancy of one H per tetrahedral site is not realistic as the H- 
H interaction in Fe is repulsive. We have calculated the H-H interaction as a function 
of distance and find strong H-H repulsion up to third-nearest neighbours. For a realistic 
representation of the statistics, up to and including third nearest-neighbour H-H spacings 
should be strongly unfavourable in our model. In the extreme case, these sites are left 
unoccupied, corresponding to infinite repulsion. The expression for the entropy derivative 
for tetrahedral interstitials in a bcc lattice with total exclusion of first nearest-neighbours, 
both first and second-nearest-neighbours, and up to third-nearest neighbours is given by 
aoureau— 'as - /vm^(^_^^^^ wnere 2; - (n-6iv)27v > dn - '^^^{n-au)^^ 



and '^^^ = -kin 

an 



n{N-3n){N-2n)^ 



respectively. In the limits of low concentration, all of the 



{N-4n)(N-5ny 

expressions for the entropy {Sq, Si, Sinn, S2nn, S^nn) are the same, the differences being of 
O(c^). Similarly, at small strains, all of the solutions for c are identical. 

A numerical example of the extent of the modification of the H-concentration in bcc Fe by 



applied hydrostatic strain is shown in Figure 13 1. In this low-concentration regime, all forms 



of the entropy are identical. In this example, the background (zero-strain) concentration 
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was taken to be equal to 3xl0~^, in line with typical densities of H in iroiP. At each 
temperature studied, /i was adjusted in order to achieve this background concentration. 
From the figure, it is clear that the H-concentration can be dramatically increased by orders 
of magnitude, even at very small strains. This increase in concentration is more pronounced 
at lower temperatures, in agreement with the low thermal desportion temperatures, ca 400 
K, indicated for H trapped at dislocation^^. Edge dislocations have a significant hydrostatic 
strain component^ so we take our results to be representative of these. With shear strain 
oriented in the plane normal to the (111) direction, we obtain the curves in Fig. [13}d. The 
effect is not as dramatic, because the solubility energy is quadratic in the shear strain [c.f. 



Eq. [22| whose coefficients change upon rotation to this coordinate system, but not the 
functional form). 
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FIG. 13: Concentration of H as a function of (a) hydrostatic strain, and (b) shear strain in the 
plane normal to the (111) direction, for a chemical potential giving a background concentration of 
3 X 10^^ H per Fe atom at each temperature. 



V. CONCLUSIONS 

We have studied the modification of elastic properties of H in bcc Fe. The modification 
quite profoundly affects the elastic moduli. Using this knowledge, we have characterised 
the strain-dependence for H-solubility in Fe. From our model we find, in agreement with 
existing experimental studieJ^, that typical strain fields found near dislocations can trap 
lots of H. We further show that the local H-concentration may be enhanced by several 
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orders of magnitude near dislocations with tensile strain components, which justifies our 
use of concentrations in the atomic % range for our calculations. In the case of tensile 
strain, we find that these areas of enriched local concentration tend to occur mainly at low 
temperatures, in the hundreds of K, which agrees with the H-trapping energies found in the 
thermal desorption results of Novak. In the case of screw dislocations, we also find trapping 
of H if the screw dislocation is oriented along the (111) direction but the concentration 
enhancement is not as great as for the edge dislocation. Screw dislocations along (111) are 
the predominant types of dislocations in bcc metals driving plasticitj^ and the concentration 
enhancement near screw dislocations does not rule out the HELP mechanism in explaining 
the mechanism of H-embrittlement of bcc Fe. 
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